1 Introduction/Abstract

This purpose of this document is to collect, analyze, and plot data associated with the mansucript titled "Neural Mechanisms of Thermosensation in the Human Parasite Strongyloides stercoralis.

The data analyzed here belong to two different classes of experiments: thermotaxis tracking assays, and thermosensory calcium imaging experiments.

1.1 Setup

2 Tax-4 and HisCl1 Experiments

These data show the results of the following experimental questions:

Role of Ss-tax-4 in negative thermotaxis behaviors by S. stercoralis iL3s
Effect of inactivating Ss-AFD on positive thermotaxis
Effect of inactivating Ss-AFD on negative thermotaxis

2.1 Data Inputs

These data were calculated by the TT_WormTracker Matlab code and are imported into RStudio from the results document generated by Matlab.

Change in temperature (\(^\circ\)C)
Worm speed (mm/s)
Distance ratio (total distance traveled/max displacement)

dat_inactivation_file <- c('~/Box Sync/Lab_Hallem/Astra/Writing/Bryant et al 20xx/Data/Tax4HisCl1_R_Import.xlsx')

if (is_empty(dat_inactivation_file)){
  dat_inactivation_file <- file.choose()
}
dat_inactivation <- read_excel(dat_inactivation_file, sheet = 1) %>%
  group_by(Experiment,ThermoMode, Phenotype) %>%
  dplyr::mutate(Experiment = factor(Experiment, levels = c("HisCl", "CRISPR"))) %>%
  dplyr::mutate(ThermoMode = factor(ThermoMode, levels = c("PT", "NT"))) %>%
  dplyr::mutate(Phenotype = factor(Phenotype, levels = c("BU", "Histamine", "NoCas9", "Ss-tax-4")))

2.2 Plot Templates

plot_tracking <-  function(dat, y, ylim, xlabel, ylabel, title = NULL, groupcolors,groupshapes){
  ggplot(data = dat, mapping = aes(x = Phenotype, y = y)) + 
    stat_boxplot(geom = "errorbar", 
                 width = 0.25,
                 size = .5) +
    geom_boxplot(width = 0.5,
                 outlier.shape = NA,
                 fatten = 1,
                 key_glyph = "rect",
                 alpha = 0.3,
                 show.legend = F,
                 aes(fill = Phenotype)) +
    geom_jitter(size = 1,
                stroke = .75,
                width = .2,
                aes(color = Phenotype, shape = Experiment)) +
    scale_y_continuous(limits = ylim, 
                       expand = expansion(add = 0)) +
    scale_color_manual(values = groupcolors,
                       aesthetics = c("colour", "fill")) +
    scale_shape_manual(values = groupshapes)  +
    xlab(xlabel) +
    ylab(ylabel) +
    ggtitle(title) +
    theme_Publication() +
    theme(aspect.ratio = 2/1,
          axis.text.x = element_text(angle= 45,
                                     hjust = 1)) 
}

2.3 AFDp::HisCl1

2.3.1 Positive Thermotaxis

p.PT.HisCl.dT<- dat_inactivation %>%
  dplyr::filter(Experiment == "HisCl") %>%
  dplyr::filter(ThermoMode == "PT") %>%
  plot_tracking(dat = ., 
                y = .$dT, 
                ylim = c(-2,6),
                xlabel = NULL,
                ylabel = "dT (C)",
                title = "PT HisCl1",
                groupcolors = c("BU" = "#5F6369", 
                                "Histamine" = "#99772C"),
                groupshapes = 21)

p.PT.HisCl.Speed<- dat_inactivation %>%
  dplyr::filter(Experiment == "HisCl") %>%
  dplyr::filter(ThermoMode == "PT") %>%
  plot_tracking(dat = ., 
                y = .$Speed, 
                ylim = c(0,0.8),
                xlabel = NULL,
                title = "PT HisCl1",
                ylabel = "Speed (mm/s)",
                groupcolors = c("BU" = "#5F6369", 
                                "Histamine" = "#99772C"),
                groupshapes = 21)

p.PT.HisCl.DR<- dat_inactivation %>%
  dplyr::filter(Experiment == "HisCl") %>%
  dplyr::filter(ThermoMode == "PT") %>%
  plot_tracking(dat = ., 
                y = .$`Distance Ratio`, 
                ylim = c(0,8),
                xlabel = NULL,
                title = "PT HisCl1",
                ylabel = "Distance Ratio",
                groupcolors = c("BU" = "#5F6369", 
                                "Histamine" = "#99772C"),
                groupshapes = 21)

ggarrange(
  p.PT.HisCl.dT,
  p.PT.HisCl.Speed,
  p.PT.HisCl.DR,
  ncol = 3,
  nrow = 1
)

2.3.2 Negative Thermotaxis

p.NT.HisCl.dT<- dat_inactivation %>%
  dplyr::filter(Experiment == "HisCl") %>%
  dplyr::filter(ThermoMode == "NT") %>%
  plot_tracking(dat = ., 
                y = .$dT, 
                ylim = c(-4,4),
                xlabel = NULL,
                ylabel = "dT (C)",
                title = "NT HisCl1",
                groupcolors = c("BU" = "#5F6369", 
                                "Histamine" = "#99772C"),
                groupshapes = 21)

p.NT.HisCl.Speed<- dat_inactivation %>%
  dplyr::filter(Experiment == "HisCl") %>%
  dplyr::filter(ThermoMode == "NT") %>%
  plot_tracking(dat = ., 
                y = .$Speed, 
                ylim = c(0,0.3001),
                xlabel = NULL,
                ylabel = "Speed (mm/s)",
                title = "NT HisCl1",
                groupcolors = c("BU" = "#5F6369", 
                                "Histamine" = "#99772C"),
                groupshapes = 21)

p.NT.HisCl.DR<- dat_inactivation %>%
  dplyr::filter(Experiment == "HisCl") %>%
  dplyr::filter(ThermoMode == "NT") %>%
  plot_tracking(dat = ., 
                y = .$`Distance Ratio`, 
                ylim = c(0,8),
                xlabel = NULL,
                ylabel = "Distance Ratio",
                title = "NT HisCl1",
                groupcolors = c("BU" = "#5F6369", 
                                "Histamine" = "#99772C"),
                groupshapes = 21)

ggarrange(
  p.NT.HisCl.dT,
  p.NT.HisCl.Speed,
  p.NT.HisCl.DR,
  ncol = 3,
  nrow = 1
)

2.4 Ss-tax-4 CRISPR

2.4.1 Negative Thermotaxis

p.NT.CRISPR.dT<- dat_inactivation %>%
  dplyr::filter(Experiment == "CRISPR") %>%
  dplyr::filter(ThermoMode == "NT") %>%
  plot_tracking(dat = ., 
                y = .$dT, 
                ylim = c(-4,4),
                xlabel = NULL,
                ylabel = "dT (C)",
                title = "NT Ss-tax-4 KO",
                groupcolors = c("NoCas9" = "#5F6369", 
                                "Ss-tax-4" = "#1A4C62"),
                groupshapes = 21)

p.NT.CRISPR.Speed<- dat_inactivation %>%
  dplyr::filter(Experiment == "CRISPR") %>%
  dplyr::filter(ThermoMode == "NT") %>%
  plot_tracking(dat = ., 
                y = .$Speed, 
                ylim = c(0,0.3001),
                xlabel = NULL,
                ylabel = "Speed (mm/s)",
                title = "NT Ss-tax-4 KO",
                groupcolors = c("NoCas9" = "#5F6369", 
                                "Ss-tax-4" = "#1A4C62"),
                groupshapes = 21)

p.NT.CRISPR.DR<- dat_inactivation %>%
  dplyr::filter(Experiment == "CRISPR") %>%
  dplyr::filter(ThermoMode == "NT") %>%
  plot_tracking(dat = ., 
                y = .$`Distance Ratio`, 
                ylim = c(0,25),
                xlabel = NULL,
                ylabel = "Distance Ratio",
                title = "NT Ss-tax-4 KO",
                groupcolors = c("NoCas9" = "#5F6369", 
                                "Ss-tax-4" = "#1A4C62"),
                groupshapes = 21)

ggarrange(
  p.NT.CRISPR.dT,
  p.NT.CRISPR.Speed,
  p.NT.CRISPR.DR,
  ncol = 3,
  nrow = 1
)

2.4.2 Group and Save Plots

ggarrange(
  p.PT.HisCl.dT,
  p.PT.HisCl.Speed,
  p.PT.HisCl.DR,
  p.NT.CRISPR.dT,
  p.NT.CRISPR.Speed,
  p.NT.CRISPR.DR,
  p.NT.HisCl.dT,
  p.NT.HisCl.Speed,
  p.NT.HisCl.DR,
  ncol = 3,
  nrow = 3
) %>%
  ggsave(filename ="HisClTax4_Plots.pdf",
         device = cairo_pdf,
         width = 10.5,
         height = 6.3,
         path = dirname(dat_inactivation_file))

3 AFD Response Properties

Thermosensory response properties of C. elegans AFD (strain = IK890), and S. stercoralis AFD (expression plasmid = pASB52).

The data reflect the responses of C. elegans AFD and S. stercoralis AFD to the following thermal stimuli:

Pseudo-fictive PT 20-34\(^\circ\)C ramp
Pseudo-fictive PTe 20-40\(^\circ\)C ramp
Pseudo-fictive PT15 12-26\(^\circ\)C ramp
Pseudo-fictive PT15e 12-32\(^\circ\)C ramp
Non-physiological PT 17-40\(^\circ\)C ramp
Non-physiological NT 20-12\(^\circ\)C ramp

3.1 Data Inputs

These data were preprocessed and calculated by TCI Matlab code.

dat_AFD_file <- c('~/Box Sync/Lab_Hallem/Astra/Writing/Bryant et al 20xx/Data/Calcium Imaging/AFDResponse_R_Import.xlsx')

if (is_empty(dat_AFD_file)){
  dat_AFD_file <- file.choose()
}
dat_AFD <- read_excel(dat_AFD_file, sheet = 1) %>%
  group_by(Species,Tc, Stimulus) %>%
  dplyr::mutate(Species = as_factor(Species)) %>%
  dplyr::mutate(Stimulus = factor(Stimulus, levels = c("pFPT", "pFPTe", "R4", "NT"))) %>%
  dplyr::mutate(Tc = factor(Tc, levels = c(23, 15)))

dat_AFD_15 <- dplyr::filter(dat_AFD, Tc == 15) %>%
  dplyr::filter(Stimulus != "R4") %>%
  dplyr::filter(Stimulus != "NT")

dat_AFD_23 <- dplyr::filter(dat_AFD, Tc == 23) %>%
  dplyr::filter(Stimulus != "R4") %>%
  dplyr::filter(Stimulus != "NT")

dat_AFD_allTc <- dat_AFD %>%
  dplyr::filter(Stimulus != "R4") %>%
  dplyr::filter(Stimulus != "NT")

dat_AFD_R4 <- dat_AFD %>%
  dplyr::filter(Stimulus == "R4")

dat_AFD_across <- dplyr::filter(dat_AFD, Tc == 23)  %>%
  dplyr::filter(Stimulus != "NT") %>%
  dplyr::mutate(Type = as_factor(case_when(
    Stimulus == "R4" ~ "Fast",
    Stimulus != "R4" ~ "Slow")))

3.2 Slow “Fictive” Stimuli

Responses of Ss-AFD and Ce-AFD to thermal stimuli based on the temperatures experienced by S. stercoralis iL3s performing postitive thermotaxis in a 20-34C temperature gradient. In these cases, ramp speeds are generally 0.025C/s.

3.2.1 Plot Templates

plot_AFD <- function(dat, y, ylim, xlabel, ylabel, stat.test, stat.y.position, title = NULL){
  ggplot(data = dat, mapping = aes(x = Species, y = y)) + 
    stat_boxplot(geom = "errorbar", 
                 width = 0.25,
                 size = .5) +
    geom_boxplot(width = 0.5,
                 outlier.shape = NA,
                 fatten = 1,
                 key_glyph = "rect",
                 alpha = 0.3,
                 show.legend = F,
                 aes(fill = Species)) +
    geom_jitter(size = 1,
                stroke = .75,
                width = .2,
                aes(color = Species, shape = Stimulus)) +
    stat_pvalue_manual(stat.test, 
                       y.position = stat.y.position,
                       tip.length = 0,
                       bracket.shorten = 0.2,
                       label = 'p.signif') + 
    scale_y_continuous(limits = ylim, 
                       expand = expansion(add = 0)) +
    scale_color_carto_d(palette = "Vivid") +
    scale_fill_carto_d(palette = "Vivid") +
    scale_shape_manual(values = c(21:23))  +
    xlab(xlabel) +
    ylab(ylabel) +
    ggtitle(title) +
    theme_Publication() +
    theme(aspect.ratio = 2/1) 
}


plot_AFD_Speciesfacet <- function(dat, y, ylim, xlabel, ylabel, title = NULL){
  ggplot(data = dat, mapping = aes(x = Tc, y = y)) + 
    stat_boxplot(geom = "errorbar", 
                 width = 0.25,
                 size = .5) +
    geom_boxplot(width = 0.5,
                 outlier.shape = NA,
                 fatten = 1,
                 key_glyph = "rect",
                 alpha = 0.3,
                 show.legend = F,
                 aes(fill = Species)) +
    geom_jitter(size = 1,
                stroke = .75,
                width = .2,
                aes(color = Species, shape = Stimulus)) +
    # stat_pvalue_manual(stat.test, 
    #                    y.position = stat.y.position, 
    #                    label = 'p.signif') + 
    facet_grid(. ~ Species)+
    scale_y_continuous(limits = ylim, 
                       expand = expansion(add = 0)) +
    scale_color_carto_d(palette = "Vivid") +
    scale_fill_carto_d(palette = "Vivid") +
    scale_shape_manual(values = c(21:23))  +
    xlab(xlabel) +
    ylab(ylabel) +
    ggtitle(title) +
    theme_Publication() +
    theme(aspect.ratio = 2/1) 
}

plot_AFD_Bins <- function(dat, y, ylim, xlabel, ylabel, title = NULL){
  ggplot(data = dat, mapping = aes(x = TempBin, y = y)) + 
    stat_boxplot(geom = "errorbar", 
                 width = 0.25,
                 size = .5) +
    geom_boxplot(width = 0.5,
                 outlier.shape = NA,
                 fatten = 1,
                 key_glyph = "rect",
                 alpha = 0.3,
                 show.legend = F,
                 aes(fill = Species)) +
    geom_jitter(size = 1,
                stroke = .75,
                width = .2,
                aes(color = Species, shape = Stimulus)) +
    # stat_pvalue_manual(stat.test, 
    #                    y.position = stat.y.position, 
    #                    label = 'p.signif') + 
    facet_grid(Species ~ Tc)+
    scale_y_continuous(limits = ylim, 
                       expand = expansion(add = 0)) +
    scale_color_carto_d(palette = "Vivid") +
    scale_fill_carto_d(palette = "Vivid") +
    scale_shape_manual(values = c(21:23))  +
    xlab(xlabel) +
    ylab(ylabel) +
    ggtitle(title) +
    theme_Publication() +
    theme(aspect.ratio = 2/1) 
}

3.2.2 Temperature Threshold

The temperature threshold for each worm. Threshold is defined as the temperature at which point the caclium response deviates from baseline (F0) response by at least 3 * STD of baseline, for the duration of time it takes for the stimulus ramp to advance 0.25C. For the slow temperature ramps (.025C/s) plotted here, this time is therefore 10 seconds.

Baseline (F0) responses are defined as the calcium responses at the followin temperatures:

For Tc = 23, F0 = 20C
For Tc = 15, F0 = 12C

Statistics are the results of Wilcoxon signed-rank tests.

3.2.2.1 Tc = 23C

xlabel <- "Species"
ylabel <- " Temperature (C)"
title <- "Threshold Temperature, Tc = 23C"
ylim <- c(10, 40)

stat.thresh.23 <- compare_means(Temperature ~ Species, 
                                data = dat_AFD_23,
                                alternative = "two.sided")

p.thresh.23 <- plot_AFD(dat_AFD_23, 
                        dat_AFD_23$Temperature, 
                        ylim,
                        xlabel, 
                        ylabel,
                        stat.thresh.23,
                        38,
                        title)

print(p.thresh.23)

3.2.2.2 Tc = 15C

xlabel <- "Species"
ylabel <- " Temperature (C)"
title <- "Threshold Temperature, Tc = 15C"
ylim <- c(10, 40)


p.thresh.allTc <- plot_AFD_Speciesfacet(dat_AFD_allTc, 
                                        dat_AFD_allTc$Temperature, 
                                        ylim,
                                        xlabel, 
                                        ylabel,
                                        title)

print(p.thresh.allTc)

# stat.thresh.15 <- compare_means(Temperature ~ Species, 
#                            data = dat_AFD_15,
#                            alternative = "two.sided")

# p.thresh.15 <- plot_AFD(dat_AFD_15, 
#          dat_AFD_15$Temperature, 
#          ylim,
#          xlabel, 
#          ylabel,
#          stat.thresh.15,
#          38,
#          title)
# 
# print(p.thresh.15)

3.2.3 Pearson Correlation

Is the calcium response above threshold linearly correlated with changes in temperature. Correlation between calcium repsonses and temperature are analyzed across the following temperature ranges:

For Tc = 23, from 25 - 34C (pF PT) or 25 - 40C (pF PTe)
For Tc = 15, from 17 - 26C (pF PT15) or 17 - 32C (pF PT15e)

3.2.3.1 Tc = 23C

ylabel <- "Pearson's R"
title <- "T vs Ca Correlation, Tc = 23C"
ylim <- c(-1.5, 1.5)
stat.pearson.23 <- compare_means(AboveThPearsonR ~ Species, 
                                 data = dat_AFD_23,
                                 alternative = "two.sided")

p.pearson.23 <- plot_AFD(dat_AFD_23, 
                         dat_AFD_23$AboveThPearsonR, 
                         ylim,
                         xlabel, 
                         ylabel,
                         stat.pearson.23,
                         1.3,
                         title)

print(p.pearson.23)

3.2.3.2 Tc = 15C

ylabel <- "Pearson's R"
title <- "T vs Ca Correlation, Tc = 15C"
ylim <- c(-1.5, 1.5)
xlabel <- "Cultivation Temp"
# stat.pearson.15 <- compare_means(AboveThPearsonR ~ Species, 
#                            data = dat_AFD_15,
#                            alternative = "two.sided")
# 
# p.pearson.15 <- plot_AFD(dat_AFD_15, 
#          dat_AFD_15$AboveThPearsonR, 
#          ylim,
#          xlabel, 
#          ylabel,
#          stat.pearson.15,
#          1.3,
#          title)
# 
# print(p.pearson.15)


p.pearson.allTc <- plot_AFD_Speciesfacet(dat_AFD_allTc, 
                                         dat_AFD_allTc$AboveThPearsonR, 
                                         ylim,
                                         xlabel, 
                                         ylabel,
                                         title)

print(p.pearson.allTc)

3.2.4 Spearman Correlation

Is the calcium response near threshold linearly correlated with changes in temperature. Correlation between calcium repsonses and temperature are analyzed across the following temperature ranges:

For Tc = 23, from 20 - 25C
For Tc = 15, from 12 - 17C

3.2.4.1 Tc = 23C

ylabel <- "Spearman's R"
title <- "T vs Ca Correlation, Tc = 23C"
ylim <- c(-1.5, 1.5)
stat.spearman.23 <- compare_means(AtThSpearmanR ~ Species, 
                                  data = dat_AFD_23,
                                  alternative = "two.sided")

p.spearman.23 <- plot_AFD(dat_AFD_23, 
                          dat_AFD_23$AtThSpearmanR, 
                          ylim,
                          xlabel, 
                          ylabel,
                          stat.spearman.23,
                          1.3,
                          title)

print(p.spearman.23)

3.2.4.2 Tc = 15C

ylabel <- "Spearman's R"
title <- "T vs Ca Correlation, Tc = 15C"
ylim <- c(-1.5, 1.5)
xlabel <- "Cultivation Temp"
# stat.pearson.15 <- compare_means(AboveThPearsonR ~ Species, 
#                            data = dat_AFD_15,
#                            alternative = "two.sided")
# 
# p.pearson.15 <- plot_AFD(dat_AFD_15, 
#          dat_AFD_15$AboveThPearsonR, 
#          ylim,
#          xlabel, 
#          ylabel,
#          stat.pearson.15,
#          1.3,
#          title)
# 
# print(p.pearson.15)


p.spearman.allTc <- plot_AFD_Speciesfacet(dat_AFD_allTc, 
                                          dat_AFD_allTc$AtThSpearmanR, 
                                          ylim,
                                          xlabel, 
                                          ylabel,
                                          title)

print(p.spearman.allTc)

3.2.5 Temperature eliciting maximal and minimal calcium response

Plots showing the temperature that corresponds to the maximum and minimal calcium response for each trace.

3.2.5.1 Max Response

Only use values from traces from the pF PT stimulus, not the extended versions.

3.2.5.1.1 Tc = 23C

ylabel <- "Temperature (C)"
title <- "Temp @ Max Response, Tc = 23C"
ylim <- c(15, 40)

temp_dat <- dplyr::filter(dat_AFD_23, Stimulus == "pFPT")

stat.max.23 <- compare_means(MaximalTemp ~ Species, 
                             data = temp_dat,
                             alternative = "two.sided")

## Adding missing grouping variables: `Tc`, `Stimulus`

p.max.23 <- plot_AFD(temp_dat, 
                     temp_dat$MaximalTemp, 
                     ylim,
                     xlabel, 
                     ylabel,
                     stat.max.23,
                     38,
                     title)

print(p.max.23)

3.2.5.1.2 Tc = 15C

ylabel <- "Temperature (C)"
title <- "Temp @ Max Response, Tc = 15C"
ylim <- c(10, 35)

temp_dat <- dplyr::filter(dat_AFD_15, Stimulus == "pFPT")

stat.max.15 <- compare_means(MaximalTemp ~ Species, 
                             data = temp_dat,
                             alternative = "two.sided")

## Adding missing grouping variables: `Tc`, `Stimulus`

p.max.15 <- plot_AFD(temp_dat, 
                     temp_dat$MaximalTemp, 
                     ylim,
                     xlabel, 
                     ylabel,
                     stat.max.15,
                     32,
                     title)

print(p.max.15)

3.2.5.2 Min Response

3.2.5.2.1 Tc = 23C

ylabel <- "Temperature (C)"
title <- "Temp @ Min Response, Tc = 23C"
ylim <- c(15, 40)

stat.min.23 <- compare_means(MinimalTemperature ~ Species, 
                             data = dat_AFD_23,
                             alternative = "two.sided")

## Adding missing grouping variables: `Tc`, `Stimulus`

p.min.23 <- plot_AFD(dat_AFD_23, 
                     dat_AFD_23$MinimalTemperature, 
                     ylim,
                     xlabel, 
                     ylabel,
                     stat.min.23,
                     38,
                     title)

print(p.min.23)

3.2.5.2.2 Tc = 15C

ylabel <- "Temperature (C)"
title <- "Temp @ Min Response, Tc = 15C"
ylim <- c(10, 35)

stat.min.15 <- compare_means(MinimalTemperature ~ Species, 
                             data = dat_AFD_15,
                             alternative = "two.sided")

## Adding missing grouping variables: `Tc`, `Stimulus`

p.min.15 <- plot_AFD(dat_AFD_15, 
                     dat_AFD_15$MinimalTemperature, 
                     ylim,
                     xlabel, 
                     ylabel,
                     stat.min.15,
                     32,
                     title)

print(p.min.15)

3.2.6 Group and Save AFD Plots

ggarrange(
  p.thresh.23,
  p.thresh.allTc,
  p.min.23,
  p.max.23,
  p.min.15,
  p.max.15,
  p.pearson.23,
  p.spearman.23,
  ncol = 2,
  nrow = 4
) %>%
  ggsave(filename ="AFD_Plots_pF.pdf",
         device = cairo_pdf,
         width = 7,
         height = 8.4,
         path = dirname(dat_AFD_file))

3.3 Fast Stimuli

Responses of Ss-AFD and Ce-AFD to thermal ramps based on those commonly used for C. elegans experiments (Colon-Ramos lab, in particular). In these cases, ramp speeds are generally 0.1C/s.

3.3.1 Plot Templates

source("theme_Publication.R")
plot_AFD <- function(dat, y, ylim, xlabel, ylabel, stat.test, stat.y.position, title = NULL){
  ggplot(data = dat, mapping = aes(x = Species, y = y)) + 
    stat_boxplot(geom = "errorbar", 
                 width = 0.25,
                 size = .5) +
    geom_boxplot(width = 0.5,
                 outlier.shape = NA,
                 fatten = 1,
                 key_glyph = "rect",
                 alpha = 0.3,
                 show.legend = F,
                 aes(fill = Species)) +
    geom_jitter(size = 1,
                stroke = .75,
                width = .2,
                aes(color = Species, shape = Stimulus)) +
    stat_pvalue_manual(stat.test, 
                       y.position = stat.y.position,
                       tip.length = 0,
                       bracket.shorten = 0.2,
                       label = 'p.signif') + 
    scale_y_continuous(limits = ylim, 
                       expand = expansion(add = 0)) +
    scale_color_carto_d(palette = "Vivid") +
    scale_fill_carto_d(palette = "Vivid") +
    scale_shape_manual(values = c(23))  +
    xlab(xlabel) +
    ylab(ylabel) +
    ggtitle(title) +
    theme_Publication() +
    theme(aspect.ratio = 2/1) 
}

plot_AFD_type <- function(dat, y, ylim, xlabel, ylabel, title = NULL){
  ggplot(data = dat, mapping = aes(x = Type, y = y)) + 
    stat_boxplot(geom = "errorbar",
                 width = 0.25,
                 size = .5) +
    geom_boxplot(width = 0.5,
                 outlier.shape = NA,
                 fatten = 1,
                 key_glyph = "rect",
                 alpha = 0.3,
                 show.legend = F,
                 aes(fill = Species)) +
    geom_jitter(size = 1,
                stroke = .75,
                width = .2,
                aes(color = Species, shape = Stimulus)) +
    # stat_pvalue_manual(stat.test, 
    #                    y.position = stat.y.position,
    #                    tip.length = 0,
    #                    bracket.shorten = 0.2,
    #                    label = 'p.signif') + 
    facet_grid(. ~ Species)+
    scale_y_continuous(limits = ylim, 
                       expand = expansion(add = 0)) +
    scale_color_carto_d(palette = "Vivid") +
    scale_fill_carto_d(palette = "Vivid") +
    scale_shape_manual(values = c(21:23))  +
    xlab(xlabel) +
    ylab(ylabel) +
    ggtitle(title) +
    theme_Publication() +
    theme(aspect.ratio = 2/1) 
}

3.3.2 Pearson Correlation

Is the calcium response linearly correlated with changes in temperature. Correlation between calcium repsonses and temperature are analyzed across the following temperature ranges:

For Tc = 23, from 25 - 40C (R4)

ylabel <- "Pearson's R"
title <- "T vs Ca Correlation, Tc = 23C, Stim = R4"
ylim <- c(-1.5, 1.5)
stat.pearson.R4 <- compare_means(AboveThPearsonR ~ Species, 
                                 data = dat_AFD_R4,
                                 alternative = "two.sided")

p.pearson.R4 <- plot_AFD(dat_AFD_R4, 
                         dat_AFD_R4$AboveThPearsonR, 
                         ylim,
                         xlabel, 
                         ylabel,
                         stat.pearson.R4,
                         1.3,
                         title)

print(p.pearson.R4)

3.3.3 Temperature Threshold

The temperature threshold for each worm. Threshold is defined as the temperature at which point the caclium response deviates from baseline (F0) response by at least 3 * STD of baseline, for the duration of time it takes for the stimulus ramp to advance 0.25C. For the fast temperature ramps (.1C/s) plotted here, this time is therefore 2.5 seconds.

Baseline (F0) responses are defined as the calcium responses at the followin temperatures:

For Tc = 23, F0 = 17C

Statistics are the results of Wilcoxon signed-rank tests.

xlabel <- "Species"
ylabel <- " Temperature (C)"
title <- "Threshold Temperature, Tc = 23C, Stim = R4"
ylim <- c(10, 40)

options(contrasts = c("contr.sum", "contr.poly"))
res.aov2<- aov(Temperature ~ Type * Species, data = dat_AFD_across)
Anova(res.aov2, type = "III")

res.aov3<-TukeyHSD(res.aov2, "Type:Species")
as_tibble(res.aov3$`Type:Species`, rownames = "comparison") %>%
  separate(col = comparison, into = c("comp1", "LS1", "comp2", "LS2")) %>%
  dplyr::filter(LS1 == LS2) %>%
  unite(temp1, comp1, LS1, sep = ":")%>%
  unite(temp2, comp2, LS2, sep = ":")%>%
  unite(comparison, temp1, temp2, sep = "-")

p.thresh.across <- plot_AFD_type(dat_AFD_across, 
                                 dat_AFD_across$Temperature, 
                                 ylim,
                                 xlabel, 
                                 ylabel,
                                 title)

print(p.thresh.across)

3.3.4 Tmax and Tmin responses

Plots showing the temperature that corresponds to the maximum and minimal calcium response for each trace. For now, just plotting Tc = 23.

3.3.4.1 Max Response

ylabel <- "Temperature (C)"
title <- "Temp @ Max Response, Tc = 23C, Stim = R4"
ylim <- c(15, 45)

stat.max.R4 <- compare_means(MaximalTemp ~ Species, 
                             data = dat_AFD_R4,
                             alternative = "two.sided")

## Adding missing grouping variables: `Tc`, `Stimulus`

p.max.R4 <- plot_AFD(dat_AFD_R4, 
                     dat_AFD_R4$MaximalTemp, 
                     ylim,
                     xlabel, 
                     ylabel,
                     stat.max.R4,
                     43,
                     title)

print(p.max.R4)

3.3.4.2 Min Response

ylabel <- "Temperature (C)"
title <- "Temp @ Min Response, Tc = 23C"
ylim <- c(15, 45)

stat.min.R4 <- compare_means(MinimalTemperature ~ Species, 
                             data = dat_AFD_R4,
                             alternative = "two.sided")

## Adding missing grouping variables: `Tc`, `Stimulus`

p.min.R4 <- plot_AFD(dat_AFD_R4, 
                     dat_AFD_R4$MinimalTemperature, 
                     ylim,
                     xlabel, 
                     ylabel,
                     stat.min.R4,
                     43,
                     title)

print(p.min.R4)

3.3.5 Is Tmin above T0?

xlabel <- "Stimulus"
ylabel <- "Proportion of Worms"
title <- "Deviation of Calcium Trace from T0"

temp_dat <- dat_AFD_across %>%
  dplyr::mutate(Dir_Tmin = factor(case_when(
    TminAboveT0 > 0 ~ "Above",
    TminAboveT0 == 0 ~ "Below"),
    levels = c("Above", "Below"))) %>%
  group_by(Species,Tc, Type, Dir_Tmin)

# Counts of various conditions
dplyr::count(temp_dat, TminAboveT0)

p.bar.TminAboveT0 <- ggplot(data = temp_dat, mapping = aes(x = Type)) + 
  geom_bar(aes(fill = Dir_Tmin), position = "fill") +
  scale_color_manual(values = c("Below" = "#5F6369", 
                                "Above" = "coral4"),
                     aesthetics = c("colour", "fill")) +
  facet_grid(. ~ Species)+
  xlab(xlabel) +
  ylab(ylabel) +
  ggtitle(title) +
  theme_Publication() +
  theme(aspect.ratio = 2/1) 

print(p.bar.TminAboveT0)

3.4 Negative Stimuli

Responses of Ss-AFD and Ce-AFD to negative thermal ramps based on those commonly used for C. elegans experiments (Colon-Ramos lab, in particular). In these cases, ramp speeds are generally -0.1C/s.

3.4.1 Pearson Correlation

Is the calcium response linearly correlated with changes in temperature. Correlation between calcium repsonses and temperature are analyzed across the following temperature ranges:

For Tc = 23, from 13 - 22C (NT)

dat_AFD_NT <- dat_AFD %>%
  dplyr::filter(Stimulus == "NT")

ylabel <- "Pearson's R"
title <- "T vs Ca Correlation, Tc = 23C, Stim = NT"
ylim <- c(-1.5, 1.5)
stat.pearson.NT <- compare_means(AboveThPearsonR ~ Species, 
                                 data = dat_AFD_NT,
                                 alternative = "two.sided")

p.pearson.NT <- plot_AFD(dat_AFD_NT, 
                         dat_AFD_NT$AboveThPearsonR, 
                         ylim,
                         xlabel, 
                         ylabel,
                         stat.pearson.NT,
                         1.3,
                         title)

print(p.pearson.NT)

3.4.2 Responses at specific temperatures

These plots show the calcium response at specific temperature bins.

3.4.2.1 Response at 13C

xlabel <- "Stimulus Type"
ylabel <- "Response Amplitude (% dR/R)"
title <- "Calcium Response at 13C"
ylim <- c(-40,40)

stat.NT.13.resp <- compare_means(ResponseSize_13C ~ Species, 
                                 data = dat_AFD_NT,
                                 alternative = "two.sided")

## Adding missing grouping variables: `Tc`, `Stimulus`

p.NT.13.resp <- plot_AFD(dat_AFD_NT, 
                         dat_AFD_NT$ResponseSize_13C, 
                         ylim,
                         xlabel, 
                         ylabel,
                         stat.NT.13.resp,
                         38,
                         title)

print(p.NT.13.resp)

3.4.3 Group and Save AFD Plots

ggarrange(
  p.thresh.across,
  p.pearson.R4,
  p.bar.TminAboveT0,
  p.min.R4,
  p.max.R4,
  p.pearson.NT,
  p.NT.13.resp,
  ncol = 2,
  nrow = 4
) %>%
  ggsave(filename ="AFD_Plots_R4+NT.pdf",
         device = cairo_pdf,
         width = 7,
         height = 8.4,
         path = dirname(dat_AFD_file))

4 Reversal Behaviors

Worm tracking and calcium imaging data for reversal behaviors exhibited by S. stercoralis iL3s.

4.1 Data Inputs

These data were preprocessed and calculated by TCI Matlab code.

dat_reversal_file <- c('~/Box Sync/Lab_Hallem/Astra/Writing/Bryant et al 20xx/Data/Reversal Behavior/Reversal_R_Import.xlsx')

if (is_empty(dat_reversal_file)){
  dat_reversal_file <- file.choose()
}
dat_reversal <- read_excel(dat_reversal_file, sheet = 1) %>%
  group_by(Gradient,Tc) %>%
  dplyr::mutate(Gradient = factor(Gradient, levels = c("High", "Low"))) %>%
  dplyr::mutate(Tc = factor(Tc, levels = c(23, 15)))

4.2 Plot Templates

plot_reversal_Range <-  function(dat, y, ylim, xlabel, ylabel, title = NULL){
  ggplot(data = dat, mapping = aes(x = Tc, y = y)) + 
    stat_boxplot(geom = "errorbar", 
                 width = 0.25,
                 size = .5) +
    geom_boxplot(width = 0.5,
                 outlier.shape = NA,
                 fatten = 1,
                 key_glyph = "rect",
                 alpha = 0.3,
                 show.legend = F,
                 aes(fill = Tc)) +
    geom_jitter(size = 1,
                stroke = .75,
                width = .2,
                aes(color = Tc, shape = Gradient)) +
    facet_grid(. ~ Gradient)+
    scale_y_continuous(limits = ylim, 
                       expand = expansion(add = 0)) +
    scale_color_carto_d(palette = "Vivid") +
    scale_fill_carto_d(palette = "Vivid") +
    scale_shape_manual(values = c(21:23))  +
    xlab(xlabel) +
    ylab(ylabel) +
    ggtitle(title) +
    theme_Publication() +
    theme(aspect.ratio = 2/1) 
}

4.3 Percent Reversals

This plot shows the percent of positive thermotaxing worms that engage in reversal behaviors, for each condition.

p.reversals.PT <- plot_reversal_Range(dat_reversal, 
                                      dat_reversal$Percent_PT_Reversing, 
                                      ylim = c(-5, 105),
                                      xlabel = "Cultivation Temperature", 
                                      ylabel = "% Reversing", 
                                      title = "% Positive Thermotaxis w/ Reversal")

print(p.reversals.PT)

This plot shows the percent of all worms that engage in reversal behaviors, for each condition.

p.reversals.All <- plot_reversal_Range(dat_reversal, 
                                       dat_reversal$Percent_All_Reversing, 
                                       ylim = c(-5, 105),
                                       xlabel = "Cultivation Temperature", 
                                       ylabel = "% Reversing", 
                                       title = "% All Worms w/ Reversal")

print(p.reversals.All)

5 AFD Adaptation

5.1 Data Inputs

dat_adapt_file <- c('~/Box Sync/Lab_Hallem/Astra/Writing/Bryant et al 20xx/Data/Calcium Imaging/Adaptation_R_Import.xlsx')

if (is_empty(dat_adapt_file)){
  dat_adapt_file <- file.choose()
}
dat_AFD_adaptation <- read_excel(dat_adapt_file, sheet = 1) %>%
  dplyr::mutate(Species = as_factor(Species)) %>%
  dplyr::mutate(Stimulus = factor(Stimulus, levels = c("Reversal", "pFPT", "pFPTe", "NT"))) %>%
  dplyr::mutate(Tc = factor(Tc, levels = c(23, 15))) %>%
  dplyr::mutate(UID = as_factor(UID)) %>%
  pivot_longer(cols = c(AdaptWindowEarly, AdaptWindowLate),
               names_to = "Bin",
               names_prefix = "AdaptWindow",
               values_to = "CaResponse") %>%
  dplyr::mutate(Bin = factor(Bin, levels = c("Early", "Late"))) %>%
  dplyr::mutate(Dir_late = factor(case_when(
    AdaptWindowLate_Cat > 0 ~ "Sig Diff",
    AdaptWindowLate_Cat < 0 ~ "Sig Diff",
    AdaptWindowLate_Cat == 0 ~ "Not Sig"),
    levels = c("Sig Diff", "Not Sig"))) %>%
  dplyr::mutate(Dir_early = factor(case_when(
    AdaptWindowEarly_Cat > 0 ~ "Sig Diff",
    AdaptWindowEarly_Cat < 0 ~ "Sig Diff",
    AdaptWindowEarly_Cat == 0 ~ "Not Sig"),
    levels = c("Sig Diff", "Not Sig"))) %>%
  group_by(Species,Tc, Stimulus, UID)

5.2 Tmax

This should reveal any adaptation in the response once temperature is holding at Tmax (for variable lengths of time).

5.2.1 Calcium Response Across T(max)

xlabel <- "TimingWindow"
ylabel <- "Median CaResponse, normalized to Tc"
title <- "Calcium Response as a function of time during T(max)"
ylim <- c(-50, 50)

p.Tmax.adapt <- dat_AFD_adaptation %>%
  dplyr::filter(Stimulus != "NT" & Species != "Ce") %>%
  ggplot(mapping = aes(x = Bin, y = CaResponse)) + 
  stat_boxplot(geom = "errorbar", 
               width = 0.25,
               size = .5) +
  geom_boxplot(width = 0.5,
               outlier.shape = NA,
               fatten = 1,
               key_glyph = "rect",
               alpha = 0.3,
               show.legend = F,
               aes(fill = Bin)) +
  geom_jitter(size = 1,
              stroke = .75,
              width = .2,
              aes(color = Bin, shape = Stimulus)) +
  facet_grid(. ~ Stimulus)+
  scale_y_continuous(limits = ylim, 
                     expand = expansion(add = 0)) +
  scale_color_carto_d(palette = "Vivid") +
  scale_fill_carto_d(palette = "Vivid") +
  scale_shape_manual(values = c(21:23))  +
  xlab(xlabel) +
  ylab(ylabel) +
  ggtitle(title) +
  theme_Publication() +
  theme(aspect.ratio = 2/1) 

print(p.Tmax.adapt)

5.2.2 Calcium Response Ss vs Ce at T(max)

xlabel <- "TimingWindow"
ylabel <- "Median CaResponse, normalized to Tc"
title <- "Ce vs Ss calcium response adaptation at Tmax = 26C"
ylim <- c(-50, 200)

p.Tmax.adapt_species <- dat_AFD_adaptation %>%
  dplyr::filter(Stimulus == "pFPT") %>%
  dplyr::mutate(Species = relevel(Species, "Ce")) %>%
  ggplot(mapping = aes(x = Bin, y = CaResponse)) + 
  stat_boxplot(geom = "errorbar", 
               width = 0.25,
               size = .5) +
  geom_boxplot(width = 0.5,
               outlier.shape = NA,
               fatten = 1,
               key_glyph = "rect",
               alpha = 0.3,
               show.legend = F,
               aes(fill = Bin)) +
  geom_jitter(size = 1,
              stroke = .75,
              width = .2,
              aes(color = Bin, shape = Stimulus)) +
  facet_grid(. ~ Species)+
  scale_y_continuous(limits = ylim, 
                     expand = expansion(add = 0)) +
  scale_color_carto_d(palette = "Vivid") +
  scale_fill_carto_d(palette = "Vivid") +
  scale_shape_manual(values = c(21:23))  +
  xlab(xlabel) +
  ylab(ylabel) +
  ggtitle(title) +
  theme_Publication() +
  theme(aspect.ratio = 2/1) 

print(p.Tmax.adapt_species)

5.2.3 Categorical

xlabel <- "Stimulus"
ylabel <- "Proportion of Worms"
title <- "Deviation of Calcium Trace from Holding Response"


p.bar.adapt.late <- dat_AFD_adaptation %>%
  dplyr::filter(Stimulus != "NT") %>%
  dplyr::filter(Species != "Ce") %>%
  ggplot(mapping = aes(x = Stimulus)) + 
  geom_bar(aes(fill = Dir_late), position = "fill") +
  scale_color_manual(values = c("Not Sig" = "#5F6369", 
                                "Sig Diff" = "coral4"),
                     aesthetics = c("colour", "fill")) +
  xlab(xlabel) +
  ylab(ylabel) +
  ggtitle(title) +
  theme_Publication() +
  theme(aspect.ratio = 2/1) 

print(p.bar.adapt.late)

5.3 Tmin

This should reveal any adaptation in the response once temperature is holding at Tmin (for variable lengths of time). This is for thermal stimuli approximating negative thermotaxis

5.3.1 Calcium Response

xlabel <- "TimingWindow"
ylabel <- "Median CaResponse, normalized to Tc"
title <- "Calcium Response as a function of time during T(min)"
ylim <- c(-50, 50)

p.Tmin.adapt <- dat_AFD_adaptation %>%
  dplyr::filter(Stimulus == "NT") %>%
  dplyr::mutate(Species = relevel(Species, "Ce")) %>%
  ggplot(mapping = aes(x = Bin, y = CaResponse)) + 
  stat_boxplot(geom = "errorbar", 
               width = 0.25,
               size = .5) +
  geom_boxplot(width = 0.5,
               outlier.shape = NA,
               fatten = 1,
               key_glyph = "rect",
               alpha = 0.3,
               show.legend = F,
               aes(fill = Bin)) +
  geom_jitter(size = 1,
              stroke = .75,
              width = .2,
              aes(color = Bin, shape = Stimulus)) +
  facet_grid(. ~ Species)+
  scale_y_continuous(limits = ylim, 
                     expand = expansion(add = 0)) +
  scale_color_carto_d(palette = "Vivid") +
  scale_fill_carto_d(palette = "Vivid") +
  scale_shape_manual(values = c(21:23))  +
  xlab(xlabel) +
  ylab(ylabel) +
  ggtitle(title) +
  theme_Publication() +
  theme(aspect.ratio = 1.3/1) 

print(p.Tmin.adapt)

Run a Mann-Whitney test of Ss-AFD (Tc = 23C) at 13C (early window) versus Ss-AFD (Tc = 23C) at 13C (late window).

temp <-  dplyr::filter(dat_AFD_adaptation, Stimulus == "NT" & Species == "Ss") 
compare_means(CaResponse ~ Bin,
              data = temp)

## Adding missing grouping variables: `Species`, `Tc`, `Stimulus`, `UID`

Compare responses of Ss-AFD (Tc = 23C) at 13C (early window) versus Ss-AFD (Tc = 15C) at 22C (early window). This is asking if there warming-mediated hyperpolarization and cooling-mediated hyperpolarizations are the same initial magnitude, with only the warming-mediated hyperpolarization showing rapid adaptation.

temp <-  dplyr::filter(dat_AFD_adaptation, Stimulus == "NT" & Species == "Ss" & Bin == "Early") %>%
  dplyr::ungroup() %>%
  dplyr::select(Stimulus,CaResponse) 

temp <- dplyr::filter(dat_AFD_adaptation, Stimulus == "Reversal" & Species == "Ss" & Bin == "Early") %>%
  dplyr::ungroup() %>%
  dplyr::select(Stimulus,CaResponse) %>%
  dplyr::bind_rows(temp) %>%
  dplyr::group_by(Stimulus)

compare_means(CaResponse ~ Stimulus,
              data = temp)

Compare responses of Ss-AFD (Tc = 23C) at 13C (late window) versus Ss-AFD (Tc = 15C) at 22C (late window). This is asking if there warming-mediated hyperpolarization and cooling-mediated hyperpolarizations are the same “adapted” magnitude. The better way to do this would be a 2-way ANoVA

temp <-  dplyr::filter(dat_AFD_adaptation, Stimulus == "NT" & Species == "Ss") %>%
  dplyr::ungroup() %>%
  dplyr::select(Stimulus,CaResponse, Bin) 

temp <- dplyr::filter(dat_AFD_adaptation, Stimulus == "Reversal" & Species == "Ss") %>%
  dplyr::ungroup() %>%
  dplyr::select(Stimulus,CaResponse, Bin) %>%
  dplyr::bind_rows(temp) 


options(contrasts = c("contr.sum", "contr.poly"))
res.aov2<- aov(CaResponse ~ Bin * Stimulus, data = temp)
Anova(res.aov2, type = "III")

res.aov3<-TukeyHSD(res.aov2, "Bin:Stimulus")
as_tibble(res.aov3$`Bin:Stimulus`, rownames = "comparison") %>%
  separate(col = comparison, into = c("comp1", "LS1", "comp2", "LS2")) %>%
  dplyr::filter(comp1 == comp2) %>%
  unite(temp1, comp1, LS1, sep = ":")%>%
  unite(temp2, comp2, LS2, sep = ":")%>%
  unite(comparison, temp1, temp2, sep = "-")

5.4 Group and Save Plots

ggarrange(
  p.reversals.PT,
  p.reversals.All,
  p.bar.adapt.late,
  p.Tmax.adapt,
  p.Tmax.adapt_species,
  p.Tmin.adapt,
  ncol = 2,
  nrow = 4
) %>%
  ggsave(filename ="Reversal+Adaptation_Plots.pdf",
         #plot = p.reversals,
         device = cairo_pdf,
         width = 7,
         height = 8.4,
         path = dirname(dat_adapt_file))

6 Ss-AFD-rGC Ecoptic Expression Analysis

Analysis and plotting for a series of experiments testing the effect of misepxresion of Ss-gcy-23.1, Ss-gcy-23.2, and Ss-gcy-23.3 in C. elegans ASE neurons, in response to the pseudo-fictive PTe 20-40\(^\circ\)C ramp.

6.1 Data Inputs

Data are preprocessed using Matlab TCI file TCI_analyzeEctopic.m, which calculates the temperature at which point an experimental trace rises above 3*rms of mean control trace for at least 6 seconds.

This calculation is run independently for each trace, and the value (in degrees C) is saved.

dat_ectopic_file <- c('~/Box Sync/Lab_Hallem/Astra/Writing/Bryant et al 20xx/Data/Calcium Imaging/Ectopic Expression/EctopicExpression_R_Import.xlsx')
#dat_file <- NULL
if (is_empty(dat_ectopic_file)){
  dat_ectopic_file <- file.choose()
}
meta_ectopic <-read_excel(dat_ectopic_file, sheet = "Metadata")
dat_ectopic <- read_excel(dat_ectopic_file, sheet = 1) %>%
  group_by(rGC,Strain)


dat_ectopic_filtered <- dat_ectopic %>%
  dplyr::filter(rGC != "ASE_ctrl") %>%
  dplyr::filter(!is.na(Thresh_temp)) %>%
  dplyr::bind_rows(dplyr::filter(dat_ectopic, rGC == "ASE_ctrl"))


dat_rGC_expression_file <- c('~/Box Sync/Lab_Hallem/Astra/Writing/Bryant et al 20xx/Data/Calcium Imaging/Ectopic Expression/SsAFDrGC_Gene_Expression_Data.xlsx')

if (is_empty(dat_rGC_expression_file)){
  dat_rGC_expression_file <- file.choose()
}
dat_rGC_expression <- read_excel(dat_rGC_expression_file, 
                                 sheet = 1,
                                 skip = 5) %>%
  dplyr::select(geneID:UniProtKB) %>%
  dplyr::select(-UniProtKB) %>%
  pivot_longer(cols = -geneID,
               names_to = c("life_stage", "sample"),
               names_sep = "-",
               values_to = "log2CPM") %>%
  dplyr::mutate(life_stage = factor(life_stage,
                                    levels = c("PF",
                                               "FLF",
                                               "iL3",
                                               "iL3a",
                                               "ppL1",
                                               "ppL3",
                                               "pfL1"))) %>%
  dplyr::group_by(geneID, life_stage)

6.2 rGC Plot Templates

source("theme_Publication.R")

plot_ectopic_noctrl <- function(dat, y, meta, ylim, xlabel,
                                ylabel, stat.test, 
                                stat.y.position,
                                title = NULL){
  ggplot(data = dat, mapping = aes(x = rGC, y = y)) + 
    
    stat_boxplot(geom = "errorbar", 
                 width = 0.25,
                 size = 0.5) +
    geom_boxplot(width = 0.5,
                 outlier.shape = NA,
                 fatten = 1,
                 key_glyph = "rect",
                 alpha = 0.3,
                 aes(fill = rGC)) +
    geom_jitter(size = 2,
                stroke = .75,
                width = .2,
                aes(color = Strain, shape = Strain)) +
    stat_pvalue_manual(stat.test, 
                       y.position = stat.y.position,
                       step.increase = 0.1,
                       tip.length = 0,
                       hide.ns = TRUE,
                       bracket.shorten = 0.2,
                       label = "p.adj.signif"
    ) + 
    scale_y_continuous(limits = ylim, 
                       expand = expansion(add = 0)) +
    xlab(xlabel) +
    ylab(ylabel) +
    ggtitle(title) +
    scale_fill_manual(values = c("darkorchid4", "cyan4", "coral4"))+
    #scale_color_carto_d(palette = "Bold") +
    scale_color_manual(values = c(rep(c("darkorchid4","cyan4", "coral4"),each = 3)))+
    scale_shape_manual(values = rep(c(21:23), 3))  +
    theme_Publication() +
    theme(aspect.ratio = 4/3,
          axis.text.x = element_text(angle= 45,
                                     
                                     hjust = 1)) 
}

6.3 Temperature Threshold

Stats are Kruskal-Wallis Test with Dunn’s Post tests for difference in response amplitude.

# ---- Temperature Threshold Plot ----

xlabel <- "Ss-AFD-rGCs"
ylabel <- " Temperature (C)"
title <- "Threshold Temperature"
ylim <- c(19, 45)

dat.subset <-dplyr::filter(dat_ectopic_filtered, rGC != "ASE_ctrl")%>%
  dplyr::mutate(Strain = as_factor(Strain)) %>%
  dplyr::mutate(rGC = as_factor(rGC))

# Kruskal-Wallis testing of response amplitude differenec
one.way <- kruskal.test(Thresh_temp ~ rGC, 
                        data = dat.subset)

# Posthoc Dunn tests with Bonferroni adjustment
post.hoc <- dunn.test::dunn.test(dat.subset$Thresh_temp, 
                                 dat.subset$rGC, 
                                 method = "bonferroni",
                                 table = FALSE,
                                 list = TRUE)

##   Kruskal-Wallis rank sum test
## 
## data: x and group
## Kruskal-Wallis chi-squared = 20.9216, df = 2, p-value = 0
## 
## 
##                            Comparison of x by group                            
##                                  (Bonferroni)                                  
## 
## List of pairwise comparisons: Z statistic (adjusted p-value)
## -------------------------------------------------------
## SSTP_0000354300 - SSTP_0000775800 :  2.900435 (0.0056)*
## SSTP_0000354300 - SSTP_0000846800 : -1.454430 (0.2187)
## SSTP_0000775800 - SSTP_0000846800 : -4.492727 (0.0000)*
## 
## alpha = 0.05
## Reject Ho if p <= alpha/2

stat.test.ectopic.thresh <- dplyr::tibble(
  comparisons = post.hoc$comparisons,
  p.adj = post.hoc$P.adjusted,
  p.adj.signif = stars.pval(post.hoc$P.adjusted)
) %>%
  mutate(p.adj.signif = case_when(p.adj.signif == " " ~ "ns",
                                  TRUE ~ as.character(p.adj.signif))) %>%
  tidyr::separate(comparisons, 
                  into = c("group1", "group2"),
                  sep = " - ")

p_ectopic_thresh<- plot_ectopic_noctrl(dat.subset, 
                                       dat.subset$Thresh_temp, 
                                       meta_ectopic,
                                       ylim,
                                       xlabel, 
                                       ylabel,
                                       stat.test.ectopic.thresh,
                                       41,
                                       title)


print(p_ectopic_thresh)

6.4 Temperature driving maximum response